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ABSTRACT 

In order to produce high dynamic range images in radio in- 
terferometry, bright extended sources need to be removed with 
minimal error. However, this is not a trivial task because the 
Fourier plane is sampled only at a finite number of points. The 
ensuing deconvolution problem has been solved in many ways, 
mainly by algorithms based on CLEAN. However, such algorithms 
that use image pixels as basis functions have inherent limitations 
and by using an orthonormal basis that span the whole image, we 
can overcome them The construction of such an orthonormal 
basis involves fine tuning of many free parameters that define the 
basis functions. The optimal basis for a given problem (or a given 
extended source) is not guaranteed. In this paper, we discuss the use 
of generalized prolate spheroidal wave functions as a basis. Given 
the geometry (or the region of interest) of an extended source and 
the sampling points on the visibility plane, we can construct the 
optimal basis to model the source. Not only does this gives us the 
minimum number of basis functions required but also the artifacts 
outside the region of interest are minimized. 

Index Terms — Radio astronomy. Radio interferometry, Decon- 
volution 

I. INTRODUCTION 

Due to the increase in computing power, high dynamic range 
imaging in radio interferometry (only limited by calibration errors) 
is achievable and is essential to produce novel scientific results. One 
of the main obstacles to reach a high dynamic range is the fact that 
only a finite number of samples of the Fourier (visibility) plane data 
is observed. Moreover, in earth rotation synthesis, these sampling 
points are not on a regular grid. If the observed field of view 
contains bright extended sources, complete removal (deconvolution) 
of such sources to reveal the faint background remains a challenge. 
The commonly used algorithm for such problems is CLEAN |[2l. 
However, as shown in 1 1] (and references therein), the usage of a set 
of image pixels (clean components) in CLEAN based algorithms 
has limitations. 

An alternative approach to using clean components is to use 
an orthonormal basis to model bright extended structure in the 
observation. Once such a model is obtained, it can be subtracted 
from the visibilities to reveal the residual (or the fainter back- 
ground). The most obvious method of deriving such a basis is to 
use the observed data itself 0. However, this has the drawback 
that the derivation of the basis can only be done once the complete 
observation is available (problematic for real time imaging etc.). 
On the other hand, we can adopt any arbitrary orthonormal basis 
to a given observation, completely independent of the data. For 
instance, in [1 1 Gauss-Hermite polynomials (shapelets) were used 
to produce high dynamic range images. The drawback in such an 
approach is the selection of free parameters (such as the number of 
basis functions, and the scale) cannot be determined in an optimal 
fashion. It requires experience of the user as well as a trial and 
error approach to fine tune such a basis for a given observation. 



Due to the noise floor, any extended source has finite support lIU, 
and the region of interest (ROI) or the support of a given extended 
source might not be optimal for a given arbitrary basis. 

In order to tackle the problem of finding the optimal basis for 
a given extended source, independent of the observed data, we 
select prolate spheroidal wave functions (PSWF) 13, ||6l. A similar 
problem has been solved in magnetic resonance imaging Q and 
we extend that result to radio interferometry in this paper. Unlike 
data derived basis functions, PSWF basis can be precomputed and 
can be reused for observations at different epochs. Furthermore, 
unlike shapelets, we suffer less from artifacts outside the ROI with 
minimal number of basis functions used. In fact, PSWF are already 
being used in radio astronomical imaging to construct a regular grid 
of sampling points in the Fourier plane |8|. We refer the reader to 
[9] for similar applications of PSWF in geoscience. 

Notation: We denote vectors in bold lowercase and matrices in 
bold uppercase. The matrix transpose, Hermitian, pseudoinverse are 
denoted by (.)^, (.)^ and (.)^ respectively. The identity matrix is 
given by I. 

II. MATHEMATICAL FOUNDATIONS 

We present the basics of interferometric imaging and the use 
of PSWF in this section. For a complete overview of radio 
interferometry, the reader is referred to ISJ. 

II-A. Interferometric imaging 
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Fig. 1. (a) Sampling points (total No) in the Fourier (visibility) 
plane, (b) Image of Nx by Ny pixels, with the support (ROI) area 
shaded. The shaded area has Nh pixels. 



We consider the visibility plane to be composed of Na sampling 
points as in Fig. [T] (a). The image has N Nx x Ny) pixels. 
However, the ROI has only Nb pixels, corresponding to the shaded 



area in Fig. [T] (b). The ROI is determined by the support of the 
source structure and the noise floor |4|. 

Let f{up,Vp) be the sampled visibihty at the p-th point in the 
visibiUty plane. Let us also denote the intensity of the ^-th pixel 



on the image as f{lq 



, These two quantities are related by the 



van Cittert-Zernike theorem HQ and can be approximated by the 
Fourier transform for images with small support as: 

Na-l 

m,^,) = E finp,v,)e'^-^'-"-+'"-"-\ qe[0,N-l] (1) 

N-1 



f{Up,Vp) = f{lq 
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p G [0, Na - 1] 
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We can represent ([T]) in vectorized form as 



f = Tf, f^T^'f 

where f = [/(/o, mo), • • • , /(/at-i , mAr-i)]^, 



(2) 



f = [/(^0,V0), . . -J{UNa-l^VNa-l)f- 

The matrix T (size Na x N) has on its ^-th row and p-th column 
^-j27r(ipuq+mpvq) ^^^^^ that^given f, an estimate for the true image 
can be formed as f = T^f but we recover the true image only 
if T^T = I. In order to satisfy this condition we need a regular 
grid of sampling points in the visibility plane. This is not the case 
in reality for radio interferometry because the sampling points are 
determined by location of individual receivers and earth rotation. 

II-B. Prolate spheroidal basis 

Our objective is to find a basis function p(/,m) (and its coun- 
terpart in the visibility plane p{u, v)) that maximizes the energy in 
the ROL The criterion for selection can be written as 



A 



Z,m)eROI 



|p(Z,m)p 



(3) 



In other words, A in ([3]) is the ratio between the energy concentrated 
in the ROI and the total energy in the image. The higher the value 
of A, we have lower sidelobes and artifacts outside the ROI. 

Let the vectorized versions of m) and p{u,v), evaluated at 
the pixels and visibility points, be p and p, respectively. Then, we 
can rewrite ^ as 



A : 



IIi.^pII 



(4) 



The selection of pixels that belong to the ROI from the full 
image vector p is done by premultiplying by I^. Thus, is an 
Nb X N matrix, constructed by removing rows (corresponding to 
pixels outside the ROI) from an x identity matrix. 

We state the main result here, which is a direct extension of JT). 
The steps required for the removal of an extended source from a 
given observation are: 

1) Construct the kernel K (size Nb by Nb) as 



K = I^T^(TT^)^Tl6 



(5) 



2) Find the eigendecomposition of K and select the eigenmodes 
with largest eigenvalues. If the i-th eigenvalue, eigenvector 
pair of K is (A^, ?7j, the i-th basis vector is 



P^ = -^T^(TT^)^TI,r7, 



(6) 



3) Represent the image as the vector b (size by 1). Decom- 
pose image b, using M basis vectors P = [po, . . . , Pm-i] 

Pm = b, m = P^b (7) 



to find the mode vecto^m (size M by 1). 
4) Find equivalent basis P = T^P in the Fourier plane and 
subtract the model from the observed data z. 



Pm 



(8) 



to get the residual data vector r. 
The complete proof of the derivation of PSWF is given in Q. For 
conipleteness, we give a sketch of proof as follows. The numerator 
of & can be written as 



lll^pf 



(a) 



III^T^pf = trace(p^TI,irT^p) (9) 



= trace(C^R"^TIbI^T^R"^C) 



Here, (a) is obtained by the substitution p = T^p and (b) is 

obtained by using ^ = R-P, where R = (TT^)^/^. It is easy to 
simplify the denominator of © by substituting lb = I in ([9]) as 



||pf = trace(C^R"^TT^R"^C) = trace(C^C) (10) 

We used the fact that R-^TT^R"^ = I to obtain (c) in 
Using 12]) and ([Tot , we can rewrite & as 

C = argmax trace(C^R"^TIbI?^T^R"^C) (H) 

C, llClNi 

The solution to ilTl is the largest eigenvalue,eigenvector pair of the 
matrix K = R"^TIbI^T^R"^ The dimension of K (Na by 
Na) makes the computation of eigendecomposition prohibitively 

expensive. In order to reduce this, we apply the transform r] — 
IJ^T^R-^C to KC = AC to get 



if T^R-'TIbT/ = Xt] 



(12) 



which gives us the kernel K in ([5]) which is of dimension Nb by 
Nb. 

II- C. Computational cost reduction 

Under the assumption Na ^ N > Nb, the cost of computing 
([Sj involves finding the pseudoinverse (TT^)^ (size Na by Na) 
whose rank is at most N. However, the rank of K will be at 
most Nb. Instead of using T, we downsample the visibility points 
to construct a matrix T of size Nd by N (Nd > N). This 
downsampling can be combined with the imaging weights used. 
Therefore, ^ and ^ are evaluated using T instead of T. However, 
when we calculate the residual in ([8]), we use the full matrix T. 

III. EXAMPLE 

We take a LOFAR (http://www.lofar.org) test observation of 
Cygnus A, at a frequency of 213 MHz as an example. The image 
of the source Cygnus A is given in Fig. Efa). The peak flux (not 
normalized) is about 20 Jy and the total is about 10 kJy. The 
objective is to subtract this source from the observed data to see the 
faint background sources. The ROI of this source, which is above 
the noise floor is given in Fig. [2 (b). The image has N — 7552 
pixels of dimension 118 by 64 and the ROI has Nb = 1826 pixels. 

The observation lasted for about 8 hours and the original 
sampling points are shown in Fig. [3] (a). The number of sampling 
points is Na ~ 15 x 10^. We downsample the sampling coverage 
(select only a subset with uniform probability) to get Nd = 8931 
which is just above N. A better approach for downsampling would 
be to combine this with imaging weights (i.e., selecting fewer short 
baselines and more long baselines). 

The dominant eigenvalues of the kernel K ^ (dimension 1826 
by 1826) is shown in Fig. IH We see that approximately the first 
100 eigenvalues are close to 1 while any eigenvalue beyond the 
200-th eigenmode is almost zero. 



In Fig. \5\ we have shown the PSWF basis vectors p corre- 
sponding to the first few largest eigenvalues. Note that the support 
of the basis vectors are entirely within the ROI, thus creating 
minimal artifacts outside the ROI. We select the first 100 basis 
functions M = 100 to construct the matrix P in (|7]). Using this, 
we decompose the image in Fig. [2] (a) to get the mode vector m. 

Using the equivalent basis in the Fourier plane, to calculate the 
residual as in®. Once the residual is obtained we make the residual 
images as shown in Fig. (6] For comparison we have also shown 
the residual image obtained using a shapelet basis function based 
deconvolution. Both methods give a residual noise level of about 1 1 
mJy far away from Cygnus A. In contrast, with traditional CLEAN 
based deconvolution, we get a residual noise of about 13 mJy. 
The shapelet based model used about 300 basis functions of many 
scales. On the other hand, the PSWF model used only about 100 
modes, which is much less. 
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Fig. 2. Cygnus A (a) Observed image at 213 MHz, deconvolved 
using CLEAN. The peak flux is about 20 Jy and to total flux is 
about 10 kJy. (b) The shaded area correspond to the ROI of 1826 
pixels. 



IV. CONCLUSIONS 

We have presented the use of prolate spheroidal wave functions 
in radio interferometric image deconvolution and have demon- 
strated its feasibility by application to a real observation. We get 
results comparable with existing techniques, but with fewer basis 
functions and with less artifacts outside the ROI. Future work will 
focus on widefield imaging and reducing the computational cost. 
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Fig. 3. The sampling points in the Fourier plane (a) Original 
coverage, for an 8 hour observation with about 15 million sampling 
points, (b) Reduced coverage by downsampling by 1500 with 
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Fig. 4. Eigenspectrum of the kernel K with the first 200 largest 
eigenvalues plotted. 



and Multichannel Signal Processing Workshop (SAM), Israel, 
pp. 69-72, 2010. 
[2] J. Hogbom, "Aperture synthesis with a non regular distribution 
of interferometer baselines," A&A Suppl, vol. 15, pp. 417-426, 
1974. 

[3] R. Levanda and A. Leshem, 'Adaptive selective sidelobe can- 
celler beamformer with applications in radio astronomy," in 
proc. IEEE 26-th Convention of Electrical and Electronics 
Engineers (lEEEI), Israel, 2010. 

[4] D. Slepian, "On bandwidth," Proc. of the IEEE, vol. 64, no. 3, 
pp. 292-300, 1976. 

[5] D. Slepian and H. O. Pollak, "Prolate spheroidal wave func- 
tions, Fourier analysis, and uncertainty-I," Bell Syst. Tech. J., 
vol. 40, no. 2, pp. 43-61, 1961. 

[6] H. J. Landau and H. O. Pollak, "Prolate spheroidal wave 
functions, Fourier analysis, and uncertainty-II," Bell Syst. Tech. 
J., vol. 40, no. 2, pp. 65-84, 1961. 

[7] M. A. Lindquist, C. Zhang, G. Glover, L. Shepp, and 




Fig. 5. Some PSWF basis vectors corresponding to the largest 
eigenvalues. The ROI is indicated by the dark curve. 




-5 -4 -3 -2 I 2 3 4 S 



(b) 

Fig. 6. Residuals after subtracting Cygnus A using (a) prolate basis 
(b) shapelet basis. Note that in (b) there is more unsubtracted flux 
towards the center of the source while in (a) the residuals are 
compact (point like). Both images have peak values of 5.5 Jy. 
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